/*
 * declareVariables.cpp
 *
 *  Created on: Jun 8, 2011
 *      Author: skyang
 */

#include "declareVariables.h"
#include <iostream>

//double locallngamma(double x, double* sgngam);
using namespace std;

void freeMemories(double * alpha, double * logGammaofAlpha, double * dev_alpha, double * dev_lnGammaofAlpha, double *deviationMatrix, double *dev_deviationMatrix){
	free(alpha);
	free(logGammaofAlpha);
	free(deviationMatrix);
	cudaFree (dev_alpha);
	cudaFree (dev_lnGammaofAlpha);
	cudaFree (dev_deviationMatrix);
	return;
}



void CudaDeclareVariablesWrapper(double * alpha, double * logGammaofAlpha, double * dev_alpha, double * dev_lnGammaofAlpha, double *deviationMatrix, double *dev_deviationMatrix, int sizeAlpha){
//	  alpha = (double *) malloc (sizeAlpha*sizeof(double));
//	  logGammaofAlpha = (double *) malloc (sizeAlpha*sizeof(double));
	  double sgngam = 0;

	  for (int i = 0; i < sizeAlpha; i++){
		  alpha[i] = 0.7 + i*0.1;
		  logGammaofAlpha[i] = locallngamma(1/alpha[i], &sgngam);
	  }

//	  cudaMalloc( (void **)&dev_alpha, sizeof(double) * sizeAlpha );
//	  cudaMalloc( (void **)&dev_lnGammaofAlpha, sizeof(double) * sizeAlpha );

	  cudaMemcpy( dev_alpha, alpha, sizeof(double) * sizeAlpha, cudaMemcpyHostToDevice );
	  cudaMemcpy( dev_lnGammaofAlpha, logGammaofAlpha, sizeof(double)*sizeAlpha, cudaMemcpyHostToDevice );

		////////////////////////////////////////////////////////////////////////////
		// Generate a deviation matrix:
		////////////////////////////////////////////////////////////////////////////

		double inc_A = 0.1;//////Iteration Increment of Alpha Parameter
		double start_A = 0.7;//////// Value to Start Alpha Iteration
		double   *DIV;
		//alpha = (double *) malloc (sizeAlpha*sizeof(double));
//		deviationMatrix = (double *) malloc(sizeAlpha*numDIV*sizeof(double));
		DIV = (double *) malloc(numDIV * sizeof(double));
		for (int i = 0; i < sizeAlpha; i++){
			  alpha[i] = start_A + i*inc_A;
		}

		for(int i = 0; i<numDIV; i++)
		{
			DIV[i] = minDIV*i;
		}

		for(int j = 0; j<numDIV; j++)
		{
			for(int i = 0; i<sizeAlpha; i++)
			{
				deviationMatrix[i + j*sizeAlpha] = pow(DIV[j],alpha[i])/windowSize;

			///	fprintf(dMat,"%lf\t",D[j][i]);
			}
			///fprintf(dMat,"\n");
		}
		///fclose(dMat);
		///////////////////////////////////////////////////////////////////////////////
		// Copy deviation matrix to GPU global memory
		///////////////////////////////////////////////////////////////////////////////

		//double *dev_deviationMatrix;
//		cudaMalloc( (void **)&dev_deviationMatrix, sizeof(double) * sizeAlpha*numDIV);
		cudaMemcpy( dev_deviationMatrix, deviationMatrix, sizeof(double) * sizeAlpha*numDIV, cudaMemcpyHostToDevice );

		/////////////////////////////////////////////////////////////////////////////////
		// Debugging Code:
		////////////////////////////////////////////////////////////////////////////////
/*		double *debug1;
		debug1= (double *) malloc (sizeAlpha*sizeof(double));
	    cudaMemcpy(debug1 , dev_alpha, sizeof(double) * sizeAlpha , cudaMemcpyDeviceToHost );

	    for (int i = 0; i < sizeAlpha; i++){
	    	cout << "alpha[" << i << "] = " << debug1[i] << endl;
	    }

	    cudaMemcpy(debug1 , dev_lnGammaofAlpha, sizeof(double) * sizeAlpha , cudaMemcpyDeviceToHost );
	    for (int i = 0; i < sizeAlpha; i++){
	    	cout << "log gamma alpha[" << i << "] = " << debug1[i] << endl;
	    }


	  free(debug1);
*/
		/////////////////////////////////////////////////////////////////////////////////
		// Finished debugging
		////////////////////////////////////////////////////////////////////////////////
	  return;


}

double locallngamma(double x, double* sgngam)
{

    double a;
    double b;
    double c;
    double p;
    double q;
    double u;
    double w;
    double z;
    int i;
    double logpi;
    double ls2pi;
    double tmp;
    double result;

    *sgngam = 0;

    *sgngam = 1;
    logpi = 1.14472988584940017414;
    ls2pi = 0.91893853320467274178;
    if( x<-34.0 )
    {
        q = -x;
        w = locallngamma(q, &tmp);
        p = floor(q);
        i = floor(p+0.5);
        if( i%2==0 )
        {
            *sgngam = -1;
        }
        else
        {
            *sgngam = 1;
        }
        z = q-p;
        if( z>0.5 )
        {
            p = p+1;
            z = p-q;
        }
        z = q*sin(ae_pi*z);
        result = logpi-log(z)-w;
        return result;
    }
    if( x<13 )
    {
        z = 1;
        p = 0;
        u = x;
        while(u>=3)
        {
            p = p-1;
            u = x+p;
            z = z*u;
        }
        while(u<2)
        {
            z = z/u;
            p = p+1;
            u = x+p;
        }
        if( z<0 )
        {
            *sgngam = -1;
            z = -z;
        }
        else
        {
            *sgngam = 1;
        }
        if( u==2)
        {
            result = log(z);
            return result;
        }
        p = p-2;
        x = x+p;
        b = -1378.25152569120859100;
        b = -38801.6315134637840924+x*b;
        b = -331612.992738871184744+x*b;
        b = -1162370.97492762307383+x*b;
        b = -1721737.00820839662146+x*b;
        b = -853555.664245765465627+x*b;
        c = 1;
        c = -351.815701436523470549+x*c;
        c = -17064.2106651881159223+x*c;
        c = -220528.590553854454839+x*c;
        c = -1139334.44367982507207+x*c;
        c = -2532523.07177582951285+x*c;
        c = -2018891.41433532773231+x*c;
        p = x*b/c;
        result = log(z)+p;
        return result;
    }
    q = (x-0.5)*log(x)-x+ls2pi;
    if(x>100000000 )
    {
        result = q;
        return result;
    }
    p = 1/(x*x);
    if( x>=1000.0 )
    {
        q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
    }
    else
    {
        a = 8.11614167470508450300*0.0001;
        a = -5.95061904284301438324*0.0001+p*a;
        a = 7.93650340457716943945*0.0001+p*a;
        a = -2.77777777730099687205*0.001+p*a;
        a = 8.33333333333331927722*0.01+p*a;
        q = q+a/x;
    }
    result = q;
    return result;

} // end locallngamma
